Probability-Changing Cluster Algorithm for Potts Models 
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We propose a new effective cluster algorithm of tuning the critical point automatically, which is an 
extended version of Swendsen-Wang algorithm. We change the probability of connecting spins of the 



same type, p = 1 — e 



-J/kBT 



in the process of the Monte Carlo spin update. Since we approach the 



canonical ensemble asymptotically, we can use the finite-size scaling analysis for physical quantities 
near the critical point. Simulating the two-dimensional Potts models to demonstrate the validity of 
the algorithm, we have obtained the critical temperatures and critical exponents which are consistent 
with the exact values; the comparison has been made with the invaded cluster algorithm. 

PACS numbers; 05.50.-f-q, 64.60.Ak, 75.fO.Hk 
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Cluster algorithms have been successfully used to 
overcome slow dynamics in the Monte Carlo simulation. 
Swendsen and Wang (SW) applied the Kasteleyn- 
Fortuin (KF) |^ representation to identify clusters of 
spins. The problem of the thermal phase transition is 
mapped onto the geometric percolation problem in the 
cluster formalism Quite recently, based on the 

cluster formalism, the multiple-percolating clusters of the 
Ising system with large aspect ratio have been studied . 

Machta et al. Q proposed another type of cluster algo- 
rithm, which is called the invaded cluster (IC) algorithm; 
this algorithm samples the critical point of a spin system 
without a priori knowledge of the critical temperature, 
ft is in contrast with the usual procedure that one makes 
simulations for various parameters to determine the crit- 
ical point. The IC algorithm has been shown to be effi- 
cient in studying various physical quantities in the critical 
region, but the ensemble is not necessarily clear. More- 
over, it has a problem of "bottlenecks" , which causes the 
broad tail in the distribution of the fraction of the ac- 
cepted satisfied bonds |Q. 

In this Letter, extending the SW algorithm we pro- 
pose a new algorithm of tuning the critical point auto- 
matically. The basic idea of our algorithm is that we 
change the probability of connecting spins of the same 
type, p ~ 1 — g-Jf^BT^ ^Yie process of the Monte 
Carlo spin update, where J is the exchange coupling. 
We decrease or increase p depending on the observation 
whether the KF clusters are percolating or not percolat- 
ing. This simple negative feed-back mechanism together 
with the finite-size scaling (FSS) property of the exis- 
tence probability (also called the crossing probability) 
Ep, the probability that the system percolates, leads to 
the determination of the critical point. Since our en- 
semble is asymptotically canonical as Ap, the amount of 
the change of p, becomes 0, the distribution functions 
of physical quantities obey the FSS; as a result, we can 
determine critical exponents using the FSS analysis. 

Let us explain the procedure for our probability- 
changing cluster (PCC) algorithm in detail. As an exam- 
ple, we consider the ferromagnetic g-state Potts model ^] 
whose Hamiltonian is given by 
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and for q — 2 this corresponds to an Ising model. The 
procedure of Monte Carlo spin update is as follows: 

1. Start from some spin configuration and some value 
of p. 

2. Construct the KF clusters using the probability p, 
and check whether the system is percolating or not. 
Update spins following the same rule as the SW al- 
gorithm, that is, flip all the spins on any KF cluster 
to one of q states. 

3. If the system is percolating (not percolating) in the 
previous test, decrease (increase) p by Ap (> 0). 

4. Back to the process 2. 

After repeating the above processes, the distribution of 
p for our Monte Carlo samples approaches the Gaussian 
distribution of which mean value is Pc{L); Pc{L) is the 
probability of connecting spins, such that the existence 
probability Ep becomes 1/2. The width of the distribu- 
tion depends on the choice of Ap in the process 2, and 
becomes in the limit of Ap 0. We should note that 
Pc{L) depends on the system size L, and Ep follows the 
FSS near the critical point, 



Ep{p,L)^X{tL^/''), t=ip,-p)/p,, 



(2) 



where pc is the critical value of p for the infinite system 
{L — > oo), and v is the correlation-length critical expo- 
nent. (As for the FSS of Ep, see Ref. for example.) 
We can estimate pc from the size dependence of Pc{L) 
using Eq. (|^) and, in turn, estimate Tc through the re- 
lation Pc = 1 ~ e"'^/^'^^". We have chosen the value of 
Ep which gives Pc{L) as 1/2 because it is the simplest. 
We may modify the update process such that this value 
is different from 1/2. 

A comment should be made here on the choice of cri- 
terion to determine percolating. Machta et al. m used 
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both the extension rule and the topological rule for their 
stopping condition. The former rule is that some cluster 
has maximum extent L in at least one of the d directions 
in rf-dimensional systems. The latter rule is that some 
cluster winds around the system in at least one of the d 
directions. We can use any rule to determine percolat- 
ing, but FSS functions for physical quantities, therefore 
Pc{L), depend on the rule. 

There is one free parameter in our algorithm; we may 
choose the difference Ap in the process 3. In the limit of 
small Ap we obtain the canonical ensemble, but it takes 
a long time to equilibrate for small Ap. Practically, we 
may start with rather large Ap, and switch to smaller Ap 
with monitoring the trail of the values of p. Small steps 
of preparation are enough for equilibration. 

In order to show the validity of the present method, we 
have made simulations for the 2c? ferromagnetic 2-state 
Potts model (Ising model) and 3-state Potts model. We 
have treated the systems with linear sizes L = 64, 128, 
256, and 512. We start with Ap = 0.01, and gradually 
decrease Ap to the final value. We have chosen this final 
value of Ap as 1/(20 x L^); the steps for preparation are 
10,000 for the largest size {L — 512). After reaching the 
final small value of Ap, we have taken 100,000 (200,000) 
Monte Carlo samples in the case of g = 2 (g = 3) with 
keeping Ap as constant. From each bond configuration, 
we have made 5 (10) spin configurations in order to get 
better statistics for magnetization in the case of g = 2 
(g = 3). We have performed several runs for each size, 
and have checked the statistical errors. 
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FIG. 1. Plot of Tc{L) (in units of J/ks) as a function 
of L"^^" for the 2d Ising model [q = 2), where u—1. The 
system sizes are L = 64, 128, 256, and 512. In the inset, the 
distribution of p, /(p), for L = 64 with the topological rule is 
shown for both the FCC and IC algorithms. Different scales 
are used for the vertical axis in the inset to express two quite 
different data. 

Let us show the results of the 2d Ising model (g = 2). 
In Fig. |l|, we plot the size dependence of Tf.{L) for both 
rules in units of J/ks- We have determined Pc{L) from 
the average of p, and calculated Tc{L) through the re- 



lation Pc(i) = 1 — e"''/'''^^^^^). In this plot, as an il- 
lustration, we have used the known value of v for the 
2d Ising model = 1). Using the least square method, 
we estimate as 1.1344(2) (1.1346(2)) for the exten- 
sion (topological) rule, which is consistent with the exact 
value, [ln(l + \/2)]^^ = 1.1346. Here, the nmnbcr in the 
parenthesis denotes the uncertainty in the last digit. We 
have used the known value of v but we may treat v as 
an unknown parameter to be determined. Assuming the 
FSS relation, Tc(L) = Tc + aL-^/", which is derived from 
Eq. (^, we may follow the three-parameter (Tc,l/z^, a) 
non-linear fitting procedure. Then, we have obtained 
(Te, l/i/) = (l. 1345(2), 1.00(4)) for the extension rule, and 
(1.1344(2), 1.04(4)) for the topological rule. Both esti- 
mates of Tc and v are consistent with the exact values. 

We plot the distribution of p, /(p), for L = 64 with 
the topological rule, as an example, in the inset of Fig. [|, 
which shows that the distribution of p is sharply peaked 
at 0.58196 with the standard deviation of 0.0004 for our 
choice of the final Ap. For comparison, we also show /(p) 
for the IC algorithm of the same size with the same rule; 
different scales are used for the vertical axis in the inset 
to express two quite different data. We notice that the 
distribution of p for the IC algorithm is far broader. A 
simple linear analysis yields that the width of the distri- 
bution for the PCC algorithm is proportional to i/Ap/a, 
where a is the value of dEp/dp at Pc{L); using the FSS 
we expect a cx L^^'^ . It should be noted that we have ob- 
tained the expected Ap- and L-dependence for the width 
of the distribution of p. We use the average value of p 
for the estimate of pc(L). Performing 10 runs, we have 
estimated Pc(64) as 0.581954 ± 0.000013; the statistical 
errors are very small. 
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FIG. 2. The energy histogram, f{E), of both the PCC and 
IC algorithms for L = 64 with the topological rule for the 2d 
Ising model {q — 2). The energy is measured in units of J. 
The energy histogram obtained by the constant-temperature 
calculation using the SW algorithm is also shown by a solid 
curve; the temperature has been chosen as Tc{L) determined 
by the PCC algorithm. 
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FIG. 3. Plot of (m^) as a function of L for the 2d Ising 
model {q = 2) in logarithmic scale. 

The resulting energy histogram, f{E), of our PCC al- 
gorithm for L — 64: with the topological rule is given in 
Fig. H The energy histogram obtained by the constant- 
temperature calculation using the SW algorithm is shown 
by a solid curve in Fig. ||; the temperature has been cho- 
sen as Tc{L) determined by the PCC algorithm. The 
energy histogram of the PCC algorithm is indistinguish- 
able from that by the constant-temperature calculation 
because of the sharp distribution of p. Thus, we may say 
that our ensemble is actually canonical for small enough 
Ap. The energy histogram of the IC algorithm, which is 
also given in Fig. |[ has broad tails for both high-energy 
and low-energy sides. Although our ensemble is almost 
canonical, there are deviations in physical quantities, in 
principle; the variance of energy, (£'^) — (i?)^, becomes 
larger with non-zero Ap, for example. We have checked 
the Ap-dependence of the systematic deviation for large 
Ap. However, the deviation of the variance of energy is 
smaller than the statistical error, 2 % for L = 64, with 
our choice of Ap. 
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FIG. 4. FSS plot of p{m) for the 2d Ising model (g=2), 
where (3/^=1/8. The rules to determine percolating are the 
extension rule (a) and the topological rule (b). 



In order to estimate another critical exponent f3, the 
magnetization exponent, we plot the average of the 
squared magnetization (m^) as a function of L in log- 
arithmic scale in Fig. Since our Monte Carlo sam- 
ples are sharply peaked ai p — Pc{L), in other words, at 
T = Tc{L), we can use the FSS relation. 
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for the estimate of (3 /v. From the slopes of the data for 
both rules, we have 13/v = 0.125(2) (0.126(2)) for the 
extension (topological) rule, which is consistent with the 
exact value, 1/8 (=0.125). 

It is quite interesting to study the distribution func- 
tion of physical quantities. We show the FSS plot of the 
distribution function p(rn) in Fig. ^, based on the FSS 
relation, 



p(HT=T.(L)~i^/V(mL'3/''). 



(4) 



The scaling plot for the extension rule (a) and that for 
the topological rule (b) are given there. The data for var- 
ious sizes are collapsed on a single curve. We have very 
good FSS behavior for both rules. Two rules give differ- 
ent tL^/^ in Eq. (||) for Ep — 1/2. It is easier to percolate 
for the extension rule compared with the topological rule. 
Therefore, Tc{L) of the extension rule is higher than that 
of the topological rule, which results in the difference in 
the FSS functions for p{m) between two rules. 
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FIG. 5. Plot of Tc{L) (in units of J/ks) as a function 
of L^^l" for the 2d 3-state Potts model, where i/=5/6. The 
system sizes are L = 64, 128, 256, and 512. 

Next turn to the case of the 3-state Potts model. The 
size dependence of T^iV) for both rules is shown in Fig. ||. 
We have used the known value of v for this plot; the ex- 
ponent V for the Id 3-state Potts model conjectured by 
the conformal field theory is 5/6. We estimate the 
extrapolated value of as 0.99490(6) (^.99494(6)) for 
the extension (topological) rule from Fig. ||. This value is 
consistent with the exact value, [ln(l-|--\/3)]^^ = 0.99497. 
The convergence is very good for the 3-state Potts model. 
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It is in contrast with the situation of the IC algorithm 
, where the convergence is not good enough maybe due 
to the problem of "bottlenecks" . 

We can estimate the critical exponent fj/v for the 3- 
state Potts model from the size dependence of {m?) as 
in the case of the Ising model. For the order parame- 
ter of the 3-state Potts model, we use the vector order 
parameter {ma:,my). The x and y components of the 
vector order parameter are obtained from the three com- 
ponents, nil, 7712 and ma, as tUx — 77ii — (l/2)(m2 -I-TO3) 
and 7ny — (V3/2)(m2 — 7713). Using the FSS rela- 
tion (0) for the 3-state Potts model, we have [3/v — 
0.131(2) (0.134(2)) for the extension (topological) rule, 
which is again consistent with the conjectured value, 2/15 
(=0.1333) |10|. We have also found nice FSS behavior for 
the order-parameter distribution functions p(m) of the 3- 
state Potts model as in the case of the Ising model. Here, 
771 stands for the absolute value of the vector order pa- 
rameter. The details will be given in a separate paper. 

Here we may remark on the computation time. We 
only need to modify the code of the SW algorithm slightly 
in order to change p. The typical computation time to 
get 10^ Monte Carlo samples for L = 64 is 471 seconds 
using the Alpha 21164A (533 MHz) machine, which is 
about 30 % longer than that for the SW algorithm, 361 
seconds. The most of time is spent on the usual pro- 
cedure of the SW algorithm, that is, the assignment of 
the cluster and the cluster spin update. With a small 
cost of computation time, we can determine the critical 
point automatically without making simulations for var- 
ious parameters. In contrast, it takes 958 seconds in the 
case of the IC algorithm to get 10^ samples for L = 64 be- 
cause one should check whether the system is percolating 
several times before getting one Monte Carlo sample. 

To summarize, we have proposed a new cluster algo- 
rithm of tuning the critical point automatically. Our al- 
gorithm is the extension of the SW algorithm |Q] , but we 
change the probability of connecting spins p in the pro- 
cess of Monte Carlo spin update. The resulting distribu- 
tion of p is sharply peaked at Pc{L). We approach the 
canonical ensemble in the limit of small Ap, which has 
been explicitly checked by the energy histogram. This 
is in contrast with the histogram of the IC algorithms 
0, which has broad tails for both high-energy and low- 
energy sides. Thus, we can use the FSS analysis for the 
physical quantities. The estimated values of the critical 
temperatures and the critical exponents for the 2d Potts 
models are consistent with the exact values. In order to 
get more accurate estimate of the critical point and crit- 
ical exponents, the FSS analysis employed by Ferrenberg 
and Landau ||l|] in a high-resolution Monte Carlo study 
is useful for the data obtained by the PCC algorithm; 
we extract u first, and with v determined quite accu- 
rately we can estimate Tc- Since the main purpose of the 
present Letter is to present a new and simple idea of the 
cluster algorithm, the refined data analysis including the 



corrections to FSS will be left to a subsequent study. 

In the present study, we have shown the application 
of the PCC algorithm to the thermal phase transition of 
Potts models, but the idea is based only on the prop- 
erty of a percolation problem. Thus, of course, we can 
use this algorithm in the study of the geometric percola- 
tion problem, that is, the q = 1 Potts model, to get the 
percolation threshold, pc, and various critical properties. 

Moreover, we can apply the PCC algorithm to any 
problem where the mapping to the cluster formalism ex- 
ists. It is straightforward to apply this algorithm to the 
diluted Ising (Potts) model. The lack of self-averaging 
has been recently discussed for the three-dimensional di- 
luted Ising model |l^Jl^. The PCC algorithm is quite 
useful for the problem with the lack of self-averaging, 
where the distribution of Tc due to the randomness is 
important, because we can tune the critical point of each 
random sample automatically. Another direction of ap- 
plication is the cluster quantum Monte Carlo simulation 
[|l^,^, and this problem is left to a future study. 
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